# -*- coding: utf-8 -*-
"""
Created on Fri Nov 25 15:51:00 2016

@author: qiyu
"""


import pandas as pd
from pandas import DataFrame
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import mode
from sklearn.cross_validation import train_test_split
from sklearn.svm import SVC
from sklearn.preprocessing import LabelEncoder
from sklearn.metrics import precision_recall_curve
import os
import pylab as pl
from sklearn.linear_model import LinearRegression 
from sklearn.linear_model import LogisticRegression
import statsmodels.api as sm
import statsmodels.formula.api as smf
import scipy.stats as stats
from scipy.stats import chisquare 
#import seaborn as sns

##导入数据
os.chdir('D:\\model\\salary')
sa = pd.read_csv('data.csv',sep='\001', encoding='utf-8', header=None)
sa.columns=['cv_id','gender','school','degree','major','location','age','work_year','industry','salary']

##数据探索
for i in sa.columns:
    if sa[i].dtypes=='object'  and i != 'cv_id':
        print sa[i].value_counts()

##填充异常值
sa['work_year'] = sa['work_year'].replace('\N',np.nan)
sa['age'] = sa['age'].replace('\N',np.nan)
sa['salary'] = sa['salary'].replace('\N',np.nan)
sa['work_year'] = sa['work_year'].astype('float')
sa['age'] = sa['age'].astype('float')
sa['salary'] = sa['salary'].astype('float')
#sa['school'] = sa['school'].astype('float')
#sa['location'] = sa['location'].astype('float')
#sa['salary'] = sa['gender'].astype('float')
#work_year异常值去除

s1=sa[sa['age']>=18] 
s2=s1[s1['age']<=40]         
s3=s2[s2['work_year']>=0]
s4=s3[s3['work_year']<=22]
sa1=s4[s4['gender']>0]

#plt.scatter(sa1.work_year,sa1.salary,alpha=0.3)
#plt.xlabel('work_year')
#plt.ylabel('salary')

     
sa1.ix[sa1.gender==u'男','gender']="1"
sa1.ix[sa1.gender==u'女','gender']="0"         

sa1.ix[sa1.degree==u'专科','degree']="1"
sa1.ix[sa1.degree==u'本科','degree']="2"
sa1.ix[sa1.degree==u'硕士','degree']="3"
sa1.ix[sa1.degree==u'博士','degree']="3" 

sa1.ix[sa1.school==1,'school']="1"
sa1.ix[sa1.school==2,'school']="1"
sa1.ix[sa1.school==3,'school']="1"
sa1.ix[sa1.school==4,'school']="2" 
sa1.ix[sa1.school==5,'school']="3"


sa1.ix[sa1.major==u'其它','major']       ="1"
sa1.ix[sa1.major==u'其他专业','major']       ="1"
sa1.ix[sa1.major==u'艺术学','major']         ="2"
sa1.ix[sa1.major==u'文化艺术','major']       ="2"
sa1.ix[sa1.major==u'工学','major']           ="3"
sa1.ix[sa1.major==u'管理学','major']         ="4"
sa1.ix[sa1.major==u'经济学','major']         ="5"
sa1.ix[sa1.major==u'法学','major']           ="6"
sa1.ix[sa1.major==u'文学','major']           ="7"
sa1.ix[sa1.major==u'理学','major']           ="8"
sa1.ix[sa1.major==u'教育学','major']         ="9"
sa1.ix[sa1.major==u'历史学','major']         ="10"
sa1.ix[sa1.major==u'土木建筑','major']       ="11"
sa1.ix[sa1.major==u'装备制造','major']       ="12"
sa1.ix[sa1.major==u'电子信息','major']       ="13"
sa1.ix[sa1.major==u'医药卫生','major']       ="14"
sa1.ix[sa1.major==u'财经商贸','major']       ="15"
sa1.ix[sa1.major==u'旅游','major']           ="16"
sa1.ix[sa1.major==u'教育与体育','major']     ="17"
sa1.ix[sa1.major==u'公共管理与服务','major'] ="18"
sa1.ix[sa1.major==u'医学','major']           ="19"
sa1.ix[sa1.major==u'农林牧渔','major']       ="20"
sa1.ix[sa1.major==u'农学','major']           ="21"
sa1.ix[sa1.major==u'军事学','major']         ="1"
sa1.ix[sa1.major==u'博士交叉学科','major']   ="1"
sa1.ix[sa1.major==u'新闻传播','major']       ="1"
sa1.ix[sa1.major==u'哲学','major']           ="1"
sa1.ix[sa1.major==u'资源环境与安全','major'] ="1"
sa1.ix[sa1.major==u'能源动力与材料','major'] ="1"
sa1.ix[sa1.major==u'交通运输','major']       ="1"
sa1.ix[sa1.major==u'公安与司法','major']     ="1"
sa1.ix[sa1.major==u'食品药品与粮食','major'] ="1"
sa1.ix[sa1.major==u'轻工纺织','major']       ="1"
sa1.ix[sa1.major==u'水利','major']           ="1"
sa1.ix[sa1.major==u'生物与化工','major']     ="1"
sa1.ix[sa1.major==u'硕士交叉学科','major']   ="1"



sa1.ix[sa1.industry==u'计算机软件','industry']   ="7 "
sa1.ix[sa1.industry==u'互联网/电子商务','industry']   ="7 "
sa1.ix[sa1.industry==u'房地产开发','industry']   ="11"
sa1.ix[sa1.industry==u'机械/设备/重工','industry']   ="3 "
sa1.ix[sa1.industry==u'电子技术/半导体/集成电路','industry']   ="7 "
sa1.ix[sa1.industry==u'汽车及零配件','industry']   ="3 "
sa1.ix[sa1.industry==u'通信/电信/网络设备','industry']   ="7 "
sa1.ix[sa1.industry==u'金融/投资/证券','industry']   ="10"
sa1.ix[sa1.industry==u'建筑/建材/工程','industry']   ="5 "
sa1.ix[sa1.industry==u'其他行业','industry']   ="20"
sa1.ix[sa1.industry==u'制药/生物工程','industry']   ="3 "
sa1.ix[sa1.industry==u'服装/纺织/皮革','industry']   ="3 "
sa1.ix[sa1.industry==u'快速消费品(食品、饮料、化妆品)','industry']   ="8 "
sa1.ix[sa1.industry==u'贸易/进出口','industry']   ="8 "
sa1.ix[sa1.industry==u'银行','industry']   ="10"
sa1.ix[sa1.industry==u'广告','industry']   ="15"
sa1.ix[sa1.industry==u'交通/运输/物流','industry']   ="6 "
sa1.ix[sa1.industry==u'批发/零售','industry']   ="8 "
sa1.ix[sa1.industry==u'仪器仪表/工业自动化','industry']   ="3 "
sa1.ix[sa1.industry==u'通信/电信运营、增值服务','industry']   ="7 "
sa1.ix[sa1.industry==u'医疗设备/器械','industry']   ="3 "
sa1.ix[sa1.industry==u'快速消费品(食品、饮料、化妆品)','industry']   ="8 "
sa1.ix[sa1.industry==u'石油/化工/矿产/地质','industry']   ="2 "
sa1.ix[sa1.industry==u'教育/培训/院校','industry']   ="16"
sa1.ix[sa1.industry==u'多元化业务集团公司','industry']   ="21"
sa1.ix[sa1.industry==u'计算机硬件','industry']   ="7 "
sa1.ix[sa1.industry==u'计算机服务（系统、数据服务，维修）','industry']   ="7 "
sa1.ix[sa1.industry==u'专业服务(咨询、人力资源、财会)','industry']   ="12"
sa1.ix[sa1.industry==u'网络游戏','industry']   ="7 "
sa1.ix[sa1.industry==u'建筑与工程','industry']   ="5 "
sa1.ix[sa1.industry==u'保险','industry']   ="10"
sa1.ix[sa1.industry==u'房地产/建筑/建材/工程','industry']   ="11"
sa1.ix[sa1.industry==u'物业管理/商业中心','industry']   ="12"
sa1.ix[sa1.industry==u'石油/化工/矿产','industry']   ="2 "
sa1.ix[sa1.industry==u'酒店/旅游','industry']   ="9 "
sa1.ix[sa1.industry==u'原材料和加工','industry']   ="3 "
sa1.ix[sa1.industry==u'餐饮业','industry']   ="9 "
sa1.ix[sa1.industry==u'会计/审计','industry']   ="12"
sa1.ix[sa1.industry==u'家居/室内设计/装潢','industry']   ="15"
sa1.ix[sa1.industry==u'法律','industry']   ="12"
sa1.ix[sa1.industry==u'专业服务（咨询，人力资源）','industry']   ="12"
sa1.ix[sa1.industry==u'医疗/护理/卫生','industry']   ="17"
sa1.ix[sa1.industry==u'环保','industry']   ="17"
sa1.ix[sa1.industry==u'中介服务','industry']   ="12"
sa1.ix[sa1.industry==u'家具/家电/玩具/礼品','industry']   ="8 "
sa1.ix[sa1.industry==u'家具/家电/工艺品/玩具/珠宝','industry']   ="8 "
sa1.ix[sa1.industry==u'影视/媒体/艺术/文化传播','industry']   ="18"
sa1.ix[sa1.industry==u'公关/市场推广/会展','industry']   ="18"
sa1.ix[sa1.industry==u'家具/家电/工艺品/玩具','industry']   ="8 "
sa1.ix[sa1.industry==u'新能源','industry']   ="14"
sa1.ix[sa1.industry==u'教育/培训','industry']   ="16"
sa1.ix[sa1.industry==u'文字媒体/出版','industry']   ="18"
sa1.ix[sa1.industry==u'互联网','industry']   ="7 "
sa1.ix[sa1.industry==u'电力/水利','industry']   ="4 "
sa1.ix[sa1.industry==u'电气/电力/水利','industry']   ="4 "
sa1.ix[sa1.industry==u'印刷/包装/造纸','industry']   ="3 "
sa1.ix[sa1.industry==u'金融','industry']   ="10"
sa1.ix[sa1.industry==u'医药/生物工程','industry']   ="3 "
sa1.ix[sa1.industry==u'快速消费品','industry']   ="8 "
sa1.ix[sa1.industry==u'外包服务','industry']   ="12"
sa1.ix[sa1.industry==u'农/林/牧/渔','industry']   ="1 "
sa1.ix[sa1.industry==u'IT服务(系统/数据/维护)','industry']   ="7 "
sa1.ix[sa1.industry==u'生活服务','industry']   ="15"
sa1.ix[sa1.industry==u'计算机服务（系统、数据服务，维修）','industry']   ="7 "
sa1.ix[sa1.industry==u'娱乐/休闲/体育','industry']   ="18"
sa1.ix[sa1.industry==u'美容/保健','industry']   ="15"
sa1.ix[sa1.industry==u'办公用品及设备','industry']   ="8 "
sa1.ix[sa1.industry==u'检测，认证','industry']   ="12"
sa1.ix[sa1.industry==u'消费品','industry']   ="8 "
sa1.ix[sa1.industry==u'航天/航空','industry']   ="19"
sa1.ix[sa1.industry==u'通信','industry']   ="7 "
sa1.ix[sa1.industry==u'建筑','industry']   ="5 "
sa1.ix[sa1.industry==u'电子技术','industry']   ="7 "
sa1.ix[sa1.industry==u'政府/公共事业','industry']   ="19"
sa1.ix[sa1.industry==u'学术/科研','industry']   ="13"
sa1.ix[sa1.industry==u'奢侈品/收藏品/工艺品/珠宝','industry']   ="8 "
sa1.ix[sa1.industry==u'批发','industry']   ="8 "
sa1.ix[sa1.industry==u'计算机硬件及网络设备','industry']   ="7 "
sa1.ix[sa1.industry==u'影视/媒体/艺术','industry']   ="18"
sa1.ix[sa1.industry==u'医疗/护理/保健/卫生','industry']   ="17"
sa1.ix[sa1.industry==u'汽车/摩托车','industry']   ="3 "
sa1.ix[sa1.industry==u'印刷/包装','industry']   ="3 "
sa1.ix[sa1.industry==u'采掘业/冶炼','industry']   ="3 "
sa1.ix[sa1.industry==u'大型设备/机电设备/重工业','industry']   ="3 "
sa1.ix[sa1.industry==u'广告/会展/公关/市场推广','industry']   ="18"
sa1.ix[sa1.industry==u'服装服饰/纺织/皮革','industry']   ="3 "
sa1.ix[sa1.industry==u'家居','industry']   ="8 "
sa1.ix[sa1.industry==u'零售/批发','industry']   ="8 "
sa1.ix[sa1.industry==u'贸易','industry']   ="8 "
sa1.ix[sa1.industry==u'制药','industry']   ="3 "
sa1.ix[sa1.industry==u'政府','industry']   ="19"
sa1.ix[sa1.industry==u'机械','industry']   ="3 "
sa1.ix[sa1.industry==u'酒店/餐饮','industry']   ="9 "
sa1.ix[sa1.industry==u'仪器仪表及工业自动化','industry']   ="3 "
sa1.ix[sa1.industry==u'信息技术和互联网(计算机软硬件、通讯)','industry']   ="7 "
sa1.ix[sa1.industry==u'服装','industry']   ="3 "
sa1.ix[sa1.industry==u'专业服务(咨询，人力资源)','industry']   ="12"
sa1.ix[sa1.industry==u'计算机服务','industry']   ="7 "
sa1.ix[sa1.industry==u'石油/石化/化工','industry']   ="2 "
sa1.ix[sa1.industry==u'计算机服务(系统、数据服务、维修)','industry']   ="7 "
sa1.ix[sa1.industry==u'基金/证券/期货/投资','industry']   ="10"
sa1.ix[sa1.industry==u'物业管理','industry']   ="15"
sa1.ix[sa1.industry==u'家居/室内设计/装饰装潢','industry']   ="15"
sa1.ix[sa1.industry==u'专业服务','industry']   ="12"
sa1.ix[sa1.industry==u'医疗设备','industry']   ="3 "
sa1.ix[sa1.industry==u'食品/饮料/烟酒/日化','industry']   ="8 "
sa1.ix[sa1.industry==u'石油','industry']   ="13"
sa1.ix[sa1.industry==u'能源/矿产/采掘/冶炼','industry']   ="2 "
sa1.ix[sa1.industry==u'政府/公共事业/非盈利机构','industry']   ="19"
sa1.ix[sa1.industry==u'非盈利机构','industry']   ="19"
sa1.ix[sa1.industry==u'奢侈品/收藏品','industry']   ="8 "
sa1.ix[sa1.industry==u'交通','industry']   ="6 "
sa1.ix[sa1.industry==u'计算机','industry']   ="7 "
sa1.ix[sa1.industry==u'跨领域经营','industry']   ="20"
sa1.ix[sa1.industry==u'通讯/电信','industry']   ="7 "
sa1.ix[sa1.industry==u'农业/渔业/林业','industry']   ="1 "
sa1.ix[sa1.industry==u'家具','industry']   ="3 "
sa1.ix[sa1.industry==u'服务业','industry']   ="15"
sa1.ix[sa1.industry==u'互联网/移动互联网/电子商务','industry']   ="7 "
sa1.ix[sa1.industry==u'酒店','industry']   ="9 "
sa1.ix[sa1.industry==u'信托/担保/拍卖/典当','industry']   ="10"
sa1.ix[sa1.industry==u'医疗','industry']   ="17"
sa1.ix[sa1.industry==u'媒体/出版/影视/文化传播','industry']   ="18"
sa1.ix[sa1.industry==u'IT服务','industry']   ="7 "
sa1.ix[sa1.industry==u'影视','industry']   ="18"
sa1.ix[sa1.industry==u'仪器仪表','industry']   ="3 "
sa1.ix[sa1.industry==u'旅游/度假','industry']   ="18"
sa1.ix[sa1.industry==u'耐用消费品(服装、纺织、家具、家电、工艺品)','industry']   ="8 "
sa1.ix[sa1.industry==u'机械制造','industry']   ="3 "
sa1.ix[sa1.industry==u'加工/制造(工业自动化、设备、零部件)','industry']   ="3 "
sa1.ix[sa1.industry==u'会计','industry']   ="12"
sa1.ix[sa1.industry==u'家具/家电','industry']   ="3 "
sa1.ix[sa1.industry==u'批发和零售','industry']   ="8 "
sa1.ix[sa1.industry==u'专业服务/咨询(财会/法律/人力资源等)','industry']   ="12"
sa1.ix[sa1.industry==u'租赁服务','industry']   ="12"
sa1.ix[sa1.industry==u'电子','industry']   ="3 "
sa1.ix[sa1.industry==u'制造(机械/设备)','industry']   ="3 "
sa1.ix[sa1.industry==u'娱乐/体育/休闲','industry']   ="18"
sa1.ix[sa1.industry==u'机械制造/机电/重工','industry']   ="3 "
sa1.ix[sa1.industry==u'公关','industry']   ="18"
sa1.ix[sa1.industry==u'文字媒体','industry']   ="18"
sa1.ix[sa1.industry==u'汽车','industry']   ="3 "
sa1.ix[sa1.industry==u'医疗/护理/美容/保健/卫生服务','industry']   ="15"
sa1.ix[sa1.industry==u'基金','industry']   ="10"
sa1.ix[sa1.industry==u'广告/会展/公关','industry']   ="18"
sa1.ix[sa1.industry==u'原材料及加工','industry']   ="3 "
sa1.ix[sa1.industry==u'仪器','industry']   ="3 "
sa1.ix[sa1.industry==u'房地产开发/建筑与工程','industry']   ="11"
sa1.ix[sa1.industry==u'印刷','industry']   ="3 "
sa1.ix[sa1.industry==u'电力','industry']   ="4 "
sa1.ix[sa1.industry==u'房地产服务','industry']   ="11"
sa1.ix[sa1.industry==u'房地产开发/建筑/建材/工程','industry']   ="11"
sa1.ix[sa1.industry==u'美容','industry']   ="15"
sa1.ix[sa1.industry==u'通讯/电信（设备，运营，增值服务）','industry']   ="7 "
sa1.ix[sa1.industry==u'娱乐','industry']   ="18"
sa1.ix[sa1.industry==u'房地产','industry']   ="11"
sa1.ix[sa1.industry==u'能源','industry']   ="14"
sa1.ix[sa1.industry==u'耐用消费品','industry']   ="8 "
sa1.ix[sa1.industry==u'航空/航天','industry']   ="19"
sa1.ix[sa1.industry==u'生物/制药/保健/医药','industry']   ="3 "
sa1.ix[sa1.industry==u'加工/制造/汽车','industry']   ="3 "
sa1.ix[sa1.industry==u'化工','industry']   ="3 "
sa1.ix[sa1.industry==u'IT服务/系统集成','industry']   ="7 "
sa1.ix[sa1.industry==u'礼品/玩具/工艺美术/收藏品/奢侈品','industry']   ="8 "
sa1.ix[sa1.industry==u'旅游','industry']   ="18"
sa1.ix[sa1.industry==u'金融/保险/证券','industry']   ="10"
sa1.ix[sa1.industry==u'检验/检测/认证','industry']   ="12"
sa1.ix[sa1.industry==u'咨询业','industry']   ="12"
sa1.ix[sa1.industry==u'仪器/仪表/工业自动化/电气','industry']   ="3 "
sa1.ix[sa1.industry==u'电子技术/半导体','industry']   ="7 "
sa1.ix[sa1.industry==u'媒体/出版','industry']   ="18"
sa1.ix[sa1.industry==u'媒体','industry']   ="18"
sa1.ix[sa1.industry==u'建筑/设计/装潢','industry']   ="5 "
sa1.ix[sa1.industry==u'化工/能源','industry']   ="2 "
sa1.ix[sa1.industry==u'通信(设备/运营/增值)','industry']   ="7 "
sa1.ix[sa1.industry==u'交通/运输','industry']   ="6 "
sa1.ix[sa1.industry==u'家电业','industry']   ="3 "
sa1.ix[sa1.industry==u'航空/航天研究与制造','industry']   ="3 "
sa1.ix[sa1.industry==u'金融(银行、风险基金)','industry']   ="10"
sa1.ix[sa1.industry==u'专业服务(咨询/财会/法律/翻译等)','industry']   ="12"
sa1.ix[sa1.industry==u'房地产开发?建筑与工程','industry']   ="11"
sa1.ix[sa1.industry==u'房地产服务(物业管理/地产经纪)','industry']   ="12"
sa1.ix[sa1.industry==u'计算机硬件/网络设备','industry']   ="7 "
sa1.ix[sa1.industry==u'快速消费品（食品/饮料/烟酒/日化）','industry']   ="8 "
sa1.ix[sa1.industry==u'IT服务（系统/数据/维护）/多领域经营','industry']   ="7 "
sa1.ix[sa1.industry==u'物流/仓储','industry']   ="6 "
sa1.ix[sa1.industry==u'房地产及中介','industry']   ="12"
sa1.ix[sa1.industry==u'旅游/酒店/餐饮服务/生活服务','industry']   ="9 "
sa1.ix[sa1.industry==u'加工制造（原料加工/模具）','industry']   ="3 "
sa1.ix[sa1.industry==u'建筑业','industry']   ="5 "
sa1.ix[sa1.industry==u'计算机/互联网/通信/电子','industry']   ="7 "
sa1.ix[sa1.industry==u'耐用消费品（服饰/纺织/皮革/家具/家电）','industry']   ="8 "
sa1.ix[sa1.industry==u'房地产/建筑','industry']   ="11"
sa1.ix[sa1.industry==u'房地产开发·建筑与工程','industry']   ="11"
sa1.ix[sa1.industry==u'行业','industry']   ="20"
sa1.ix[sa1.industry==u'不详','industry']   ="20"
sa1.ix[sa1.industry==u'计算机服务(系统、数据服务，维修)','industry']   ="7 "

for i in sa1.columns:
    if sa1[i].dtypes=='object'  and i != 'cv_id':
        print sa1[i].value_counts()

sa1['gender'] = sa1['gender'].astype('float')
sa1['degree'] = sa1['degree'].astype('float')
sa1['school'] = sa1['school'].astype('float')
sa1['location'] = sa1['location'].astype('float')
sa1['work_year'] = sa1['work_year'].astype('float')
sa1['salary'] = sa1['salary'].astype('float')
sa1['major'] = sa1['major'].astype('float')
sa1['industry'] = sa1['industry'].astype('float')    
    
#看整体分布
sa1.describe() 
# 查看每一列的标准差
print sa1.std()

#相关性分析：
d1=sa1.corr() #相关系数矩阵，即给出了任意二个变量之间的相关系数

plt.hist(sa1['gender'])
plt.show()
plt.hist(sa1['age'])
plt.show()
plt.hist(sa1['degree'])
plt.show()
plt.hist(sa1['school'])
plt.show()
plt.hist(sa1['location'])
plt.show()
pl.hist(sa1['work_year'])
pl.show()
plt.hist(sa1['major'])
plt.show()
plt.hist(sa1['salary'])
plt.show()

# 频率表，表示work_year与salary的值相应的数量关系
print pd.crosstab(sa1['salary'], sa1['gender'], rownames=['salary'])
print pd.crosstab(sa1['salary'], sa1['school'], rownames=['salary'])
print pd.crosstab(sa1['salary'], sa1['degree'], rownames=['salary'])
print pd.crosstab(sa1['salary'], sa1['location'], rownames=['salary'])
print pd.crosstab(sa1['salary'], sa1['work_year'], rownames=['salary'])
a=pd.crosstab(sa1['salary'], sa1['work_year'], rownames=['salary'])
a1
## plot all of the column
sa1.hist()
pl.show()


#设虚拟变量
dummy_ranks1 = pd.get_dummies(sa1['gender'], prefix='gender')
print dummy_ranks1.head()
cols_to_keep = ['salary','work_year']
x1 = sa1[cols_to_keep].join(dummy_ranks1.ix[:,'gender_0.0':])

dummy_ranks2 = pd.get_dummies(sa1['school'], prefix='school')
print dummy_ranks2.head()
cols_to_keep1 = ['salary','work_year','gender_0.0','gender_1.0']
x2 = x1[cols_to_keep1].join(dummy_ranks2.ix[:,'school_1.0':])

dummy_ranks3 = pd.get_dummies(sa1['degree'], prefix='degree')
print dummy_ranks3.head()
cols_to_keep2 = ['salary','work_year','gender_0.0','gender_1.0','school_1.0','school_2.0','school_3.0']
x3 = x2[cols_to_keep2].join(dummy_ranks3.ix[:,'degree_1.0':])

dummy_ranks4 = pd.get_dummies(sa1['location'], prefix='location')
print dummy_ranks4.head()
cols_to_keep3 = ['salary','work_year','gender_0.0','gender_1.0','school_1.0','school_2.0','school_3.0','degree_1.0','degree_2.0','degree_3.0']
x4 = x3[cols_to_keep3].join(dummy_ranks4.ix[:,'location_1.0':])

dummy_ranks5 = pd.get_dummies(sa1['major'], prefix='major')
print dummy_ranks5.head()
cols_to_keep4 = ['salary','work_year','gender_0.0','gender_1.0','school_1.0','school_2.0','school_3.0','degree_1.0','degree_2.0','degree_3.0','location_1.0','location_2.0','location_3.0','location_4.0']
x5 = x4[cols_to_keep4].join(dummy_ranks5.ix[:,'major_1.0':])

dummy_ranks6 = pd.get_dummies(sa1['industry'], prefix='industry')
print dummy_ranks6.head()
cols_to_keep5 = ['salary','work_year','gender_0.0','gender_1.0','school_1.0','school_2.0','school_3.0','degree_1.0','degree_2.0','degree_3.0','location_1.0','location_2.0','location_3.0','location_4.0','major_1.0','major_2.0','major_3.0','major_4.0','major_5.0','major_6.0','major_7.0','major_8.0','major_9.0','major_10.0','major_11.0','major_12.0','major_13.0','major_14.0','major_15.0','major_16.0','major_17.0','major_18.0','major_19.0','major_20.0','major_21.0']  
x6 = x5[cols_to_keep5].join(dummy_ranks6.ix[:,'industry_1.0':])

#cols_to_keep4 = ['work_year*industry','gender_1.0','gender_2.0','school_2.0','school_3.0','location_2.0','location_3.0','location_4.0','major_2.0','major_3.0','major_4.0','major_5.0','major_6.0','major_7.0','major_8.0','major_9.0','major_10.0','major_11.0','major_12.0','major_13.0','major_14.0','major_15.0','major_16.0','major_17.0','major_18.0','major_19.0','major_20.0','major_21.0']  
#col=['work_year*industry']
#x11=sa1[col]

#est1=smf.ols(formula='salary ~ gender,school,work_year*industry', data=sa1).fit()
#est1.summary()

x6=sm.add_constant(x6)
y = sa1['salary'] 
y = sa1.salary

   

#x_train,x_test,y_train,y_test = train_test_split(x, y,test_size=0.3, random_state=1)
x_train,x_test = train_test_split(x6,test_size=0.3, random_state=1)
d1=(x_train, y_train)
#x_train,x_test,y_train,y_test = train_test_split(X, y,test_size=0.3, random_state=1)

#lr
regr = LinearRegression()
model=regr.fit(x_train, y_train)  
y_pred1 = regr.predict(x_test)
#
logr = LogisticRegression()
model=logr.fit(x_train, y_train) 

#最小二乘法
#sa1["intercept"]=1.0
est=sm.OLS(y_train,x_train).fit()  
est.summary()
y_pred = est.predict(x_test)


est=smf.ols(formula='salary ~g1*i7+s1*i7+s2*i7+l1*d2+l1*d3+l2*d2+l1*d3+l3*d2+l3*d3+l1*i7+l2*i7+l3*i7+l1*i10+l2*i10+l3*i10+l1*i11+l2*i11+l3*i11+work_year+s1+s2+d2+d3+m1+m2+m3+m4+m5+m6+m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18+m19+m20+i1+i2+i3+i4+i5+i6+i8+i9+i12+i13+i14+i15+i16+i17+i18+i19+i21', data=x_train).fit()
est.summary()

y_pred = est.predict(x_test)
print type(y_pred),type(x_test)  
print len(y_pred),len(x_test)  
print y_pred.shape,x_test.shape  
from sklearn import metrics  
import numpy as np  
sum_mean=0  
for i in range(len(y_pred)):  
    sum_mean+=(y_pred[i]-x_test.values[i])**2  
sum_erro=np.sqrt(sum_mean/151252)  
# calculate RMSE by hand  
print "RMSE by hand:",sum_erro 



#多元回归模型
est1=smf.OLS(y_train,x_train).fit()
est1.summary()
y_pred1 = est1.predict(x_test)

#看看训练模型都有那些属性和方法
print dir(est)

#均方根误差
print type(y_pred),type(y_test)  
print len(y_pred),len(y_test)  
print y_pred.shape,y_test.shape  
from sklearn import metrics  
import numpy as np  
sum_mean=0  
for i in range(len(y_pred)):  
    sum_mean+=(y_pred[i]-y_test.values[i])**2  
sum_erro=np.sqrt(sum_mean/151252)  
# calculate RMSE by hand  
print "RMSE by hand:",sum_erro 
5.8970841375


print type(y_pred1),type(y_test)  
print len(y_pred),len(y_test)  
print y_pred.shape,y_test.shape  
from sklearn import metrics  
import numpy as np  
sum_mean=0  
for i in range(len(y_pred1)):  
    sum_mean+=(y_pred1[i]-y_test.values[i])**2  
sum_erro=np.sqrt(sum_mean/151252)  
# calculate RMSE by hand  
print "RMSE by hand:",sum_erro
5.89885740207
  
print y_pred1
print regr.intercept_  
print regr.coef_  

x_test:
print regr.intercept_ 
5.76397793398

print regr.coef_ 
[ 2.21919124  1.36637632 -0.3982197  -0.64193637  0.32427349 -0.02681687
 -0.00661216]

#zip(feature_cols, regr.coef_) 
#zip(feature_cols, linreg.coef_) 
error=regr.predict(x)-y
model.score(x_test,y_test)
0.20523422864893293

#y_pred2 = regr.predict(x) 
#print y_pred2
#print regr.intercept_  
#print regr.coef_  



#均方根误差
print type(y_pred1),type(y_test)  
print len(y_pred1),len(y_test)  
print y_pred1.shape,y_test.shape  
from sklearn import metrics  
import numpy as np  
sum_mean=0  
for i in range(len(y_pred1)):  
    sum_mean+=(y_pred1[i]-y_test.values[i])**2  
sum_erro=np.sqrt(sum_mean/151372)  
# calculate RMSE by hand  
print "RMSE by hand:",sum_erro 


#把major作为虚拟变量
-2.40579044275
[  2.12924600e+00   2.04214100e-01   2.27568999e+00  -3.06680513e-01
  -6.98248617e-01   2.01909653e-01   6.80681695e-01  -5.18076845e-02
  -4.34346212e-01  -1.48398486e-01  -3.53048504e-01  -1.16298559e-01
  -2.38180863e-01  -4.48969684e-01  -6.79544948e-01  -1.38857240e-03
  -8.46260304e-01   8.37653032e-01   4.62026006e-01  -1.12536633e-01
  -1.45710309e-01  -2.65412316e-01  -5.55918845e-01   2.54553818e-01
   2.93642688e+00  -3.11431478e-01  -3.57143337e-01]
#RMSE
5.78037391594 
#R方
0.20992064811553646

#把industry作为虚拟变量
-0.43815694654
[ 2.0482266   0.20257253  2.20060706 -0.28883481 -0.71453762  0.20953425
  0.45945512 -0.01176109 -0.42242216 -0.25695993 -0.42864479 -0.02549672
 -0.21800556 -0.52083225 -0.66202835 -0.72024281 -0.65298859  0.77973813
  0.64323182 -0.01021809  0.07552004 -0.06004469 -0.51039215  0.73168619
  2.64953318 -0.61314758 -0.25762294 -2.78331986 -2.43516117 -3.00071365
 -2.96113737 -3.82449023 -1.63738585 -2.77137941 -3.0829568  -0.34611681
  0.35332123 -2.09352022 -4.49272369 -0.91735203 -1.47548843 -3.25581177
 -2.62301154 -2.29161641 -3.80679494 -1.33074167]

#RMSE
5.96607531914
#R方
0.23229746571657073


#r, p=stats.pearsonr(x,y) 
#chisq, p = chisquare(y_test, y_pred1)  
#print chisq  
#print p  